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Unstructured  LES  of  reacting  multiphase  flows  in 
realistic  gas  turbine  combustors 

By  Frank  Ham,  Sourabh  Apte,  Gianluca  Iaccarino,  Xiaohua  Wu,  Marcus 
Herrmann,  George  Constantinescuf,  Krishnan  Maheshij:  and  Parviz  Moin 


1.  Motivation  and  objectives 

As  part  of  the  Accelerated  Strategic  Computing  Initiative  (ASCI)  program,  an  accu- 
rate and  robust  simulation  tool  is  being  developed  to  perform  high-fidelity  LES  studies 
of  multiphase,  multiscale  turbulent  reacting  flows  in  aircraft  gas  turbine  combustor  con- 
figurations using  hybrid  unstructured  grids.  In  the  combustor,  pressurized  gas  from  the 
upstream  compressor  is  reacted  with  atomized  liquid  fuel  to  produce  the  combustion 
products  that  drive  the  downstream  turbine.  The  Large  Eddy  Simulation  (LES)  ap- 
proach is  used  to  simulate  the  combustor  because  of  its  demonstrated  superiority  over 
RANS  in  predicting  turbulent  mixing,  which  is  central  to  combustion. 

CDP  is  the  flagship  LES  code  being  developed  by  the  combustor  group  to  perform  LES 
of  reacting  multiphase  flow  in  complex  geometry.  CDP  is  named  after  the  late  Charles 
David  Pierce  (1969-2002)  who  made  several  lasting  contributions  to  the  LES  of  react- 
ing flows  and  to  this  ASCI  program.  CDP  is  a parallel  unstructured  finite-volume  code 
developed  specifically  for  LES  of  variable  density  low  Mach-number  flows.  It  is  written 
in  Fortran  90,  uses  MPI,  and  has  integrated  combustion  and  spray  modules.  In  the  first 
five  years  of  ASCI  (1997-2002),  the  principle  focus  of  the  combustor  group  was  the  de- 
velopment and  validation  of  an  entirely  new  numerical  method  with  the  characteristics 
necessary  for  simultaneously  accurate  and  robust  LES  on  unstructured  grids.  These  com- 
peting ends  were  both  achieved  by  developing  a method  around  the  principle  of  discrete 
kinetic  energy  conservation  (Mahesh  et  al.  2003).  This  numerical  algorithm  has  now  been 
extended  to  variable  density,  low-Mach  number  multiphase  reacting  flows. 

In  March  of  2003,  with  the  underlying  numerics  established  and  validated  in  a variety 
of  flows,  a major  redesign  and  rewrite  of  the  code  was  initiated.  Called  CDP-a,  this 
new  version  of  the  code  will  have  all  the  capabilities  of  the  current  code  along  with 
considerable  improvements,  all  of  which  are  considered  critical  to  achieving  the  stated 
ASCI  goal  of  a “major  advance  in  engine  simulation  technology.” 

This  paper  summarizes  the  accomplishments  of  the  combustor  group  over  the  past 
year,  concentrating  mainly  on  the  two  major  milestones  achieved  this  year: 

• Large  scale  simulation:  A major  rewrite  and  redesign  of  the  flagship  unstructured 
LES  code  has  allowed  the  group  to  perform  large  eddy  simulations  of  the  complete 
combustor  geometry  (all  18  injectors)  with  over  100  million  control  volumes. 

• Multi-physics  simulation  in  complex  geometry:  The  first  multi-physics  simulations 
including  fuel  spray  breakup,  coalescence,  evaporation,  and  combustion  are  now  being 
performed  in  a single  periodic  sector  (1  /18t/l)  of  an  actual  Pratt  & Whitney  combustor 
geometry. 
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2.  New  CDP  Code  Development:  CDP-a 

In  March  of  2003,  development  of  a new  version  of  CDP,  CDP-a,  was  initiated  with 
the  goal  of  performing  large-scale  simulation  of  the  full  Pratt  & Whitney  combustor. 
When  completed,  CDP=a  will  have  all  the  capabilities  of  the  existing  CDP  code,  with 
the  following  improvements: 

• Parallel  preprocessing  to  allow  unstructured  simulations  of  100  million  control  vol- 
umes or  more, 

• Faster,  more  scalable  solvers  based  on  geometric  and/or  algebraic  multigrid  meth- 
ods, 

• Particle/mesh  load  balancing  capabilities, 

• Integrated  postprocessing  capabilities,  including  the  parallel  computation  and  writ- 
ing of  plane  and  isosurface  cuts  of  instantaneous  or  statistical  data,  and  the  seeding  and 
tracking  of  massless  Lagrangian  particles, 

• Generally  improved  code  modularity  and  design,  allowing  users  several  levels  of 
interaction  with  the  code,  and  rapid  implementation  and  testing  of  new  models  and 
capabilities, 

• Regression  testing  and  version  control. 

At  the  time  of  writing,  the  parallel  preprocessor  has  been  completed,  and  a cold  flow 
(incompressible)  version  of  CDP-o  has  been  completed  and  verified  against  the  old  code 
in  a variety  of  incompressible  flow  calculations.  A Fortran  90  interface  to  the  algebraic 
multigrid  solver  of  LLNL’s  CASC  group  (Falgout  et  al.  2002)  has  been  added  as  an  option 
for  solving  the  Poisson  system,  resulting  in  a significant  overall  speedup,  although  at  a 
substantial  memory  overhead.  Particle/mesh  load  balancing  capabilities  have  been  added 
and  tested  using  model  particle  distributions.  Details  of  these  capabilities  are  described 
in  the  following  subsections. 

2.1.  Parallel  Preprocessing 

CDP  requires  as  input  Np  separate  grid  partition  files,  where  Np  is  the  number  of  pro- 
cessors that  will  be  used  for  the  LES.  These  partition  files  describe  the  grid  (nodes,  faces, 
control  volumes)  associated  with  a particular  partition,  and  also  contain  the  connectivity 
graph.  The  preparation  of  these  Np  partition  files  is  called  preprocessing,  and  consists  of 
four  distinct  steps: 

• geometry  definition 

• grid  generation 

• grid  partitioning 

• reordering/redistribution 

Up  to  now,  the  combustor  group  simulations  have  used  a single-processor  preprocessing 
strategy  to  generate  and  preprocess  hybrid  unstructured  grids  up  to  a maximum  size 
of  about  14  M control  volumes  (cv’s).  This  is  approximately  the  maximum  size  that 
can  fit  within  the  8 GByte  physical  memory  of  the  high-end  workstation  used  for  grid 
generation.  At  this  size,  the  grid  is  extremely  cumbersome  to  inspect  for  quality,  and  a 
more  reasonable  limit  is  about  5 M cv’s. 

To  preprocess  grids  of  order  100  M cv’s,  it  was  necessary  to  develop  a new  parallel 
preprocessor.  Some  recent  aircraft  simulations  in  the  literature  have  used  unstructured 
tetrahedral  grids  close  to  this  size  (25  to  60  M vertices)  (Mavriplis  & Pirzadeh  1999). 
These  grids  were  generated  using  a two  step  approach.  First  an  unstructured  grid  of 
several  million  vertices  was  generated  and  partitioned  on  a work-station.  These  parti- 
tions were  subsequently  refined  using  homothetic  refinement  on  a supercomputer.  This 
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Preprocessing 
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DONE  DONE 


Figure  1.  Flowchart  for  the  Parallel  Preprocessor  based  on  Fluent’s  GAMBIT  software  for 
mesh  generation,  and  our  own  preprocessing  program  called  PPAREM  (Parallel  PArtition  and 
REorder  Mesh). 


approach  significantly  reduces  the  cost  of  the  initial  grid  generation,  but  gives  up  in- 
dividual element  control.  This  approach  also  requires  a close  coupling  to  the  problem 
geometry  during  the  refinement  of  elements  along  the  boundary  to  ensure  the  grid  re- 
mains boundary-conforming. 

The  present  combustor  simulations  contain  substantially  more  geometric  complexity 
than  the  aircraft  wing-body  geometries  typically  used  for  simulation.  Consequently,  en- 
suring boundary  conformity  during  refinement  presents  a significant  challenge.  Addition- 
ally, there  are  many  regions  of  the  geometry  where  control  of  the  final  grid  distribution 
is  desired,  such  as  the  geometric  details  of  the  injectors  or  dilution  holes.  Homothetic 
refinement  without  subsequent  grid  smoothing  can  introduce  step-changes  in  the  grid 
spacing  that  may  adversely  affect  LES  accuracy.  For  these  reasons,  the  present  parallel 
preprocessor  was  developed  based  on  the  idea  of  decomposing  the  combustor  geometry 
into  several  geometrically-separated  regions,  and  then  separately  generating  the  desired 
grid  in  each  of  these  regions,  ensuring  only  that  shared  surfaces  have  identical  meshes. 
These  meshes  are  then  reassembled  on  the  parallel  computer  using  fast  octree  searches 
on  the  common  face  geometry. 

Figure  1 illustrates  the  process  of  mesh  generation  using  the  CDP  parallel  preproces- 
sor. Grid  generation  is  still  performed  on  a single  high-end  workstation  using  Fluent’s 
commercial  mesh  generation  software,  GAMBIT.  For  mesh  sizes  less  than  about  5 mil- 
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Table  1 . Parallel  preprocessing  of  35  M cv  grid 


exact  cv  count 

35,140,896 

number  of  msh  files,  Nm 

18 

number  of  final  partitions,  Np 

352 

I/O: 

input  msh  files 

3.72  GB 

output  partition  files 

3.01  GB 

Wall  clock  time: 

PPAREM  on  Nm  processors: 

read/parse  msh  files 

504  s 

octree  reconnect  geometry 

28  s 

Np  partition  (ParMETIS) 

45  s 

write  Np  partition  files 

1436  s 

Total  wall  clock  time 

2013  s (34  min) 

lion  control  volumes,  the  preprocessing  is  done  in  an  entirely  serial  way,  much  as  before. 
For  larger  meshes,  however,  the  geometric  decomposition  approach  described  earlier  is 
used,  and  GAMBIT  is  used  to  generate  a number  of  separate  mesh  files.  Our  parallel 
preprocessing  program  PPAREM  (Parallel  PArtition  and  REorder  Mesh)  then  reads  and 
reconnects  the  global  mesh,  produces  a high-quality  partition  that  is  independent  of  the 
present  geometric  decomposition  using  ParMETIS,  and  reorders  and  cooperatively  writes 
the  partition  files  using  the  MPI-I/O  routines.  Cooperative  writing  of  the  partition  files 
represented  a challenging  programming  task,  but  is  the  only  theoretically  scalable  way 
to  write  the  final  partition. 

The  parallel  preprocessor  was  used  to  prepare  two  large  grids  on  LLNL’s  Frost  for  the 
cold  flow  full  combustor  simulations  described  in  the  next  subsection.  Grid  sizes  were 
approximately  35  M and  100  M cv’s  respectively.  Due  to  the  memory  requirement  of 
ParMETIS  (approx.  0.4  GB  per  million  cv’s  out  of  the  total  requirement  of  0.5  GB  per 
million  CV’s  for  PPAREM),  it  was  necessary  in  the  case  of  the  100  M cv  grid  to  produce 
an  intermediate  partition  using  a low-memory  but  poor  quality  partitioning  algorithm 
(simply  binning  the  cv’s  based  on  global  index).  PPAREM  was  then  run  a second  time 
on  a larger  number  of  processors  using  this  poor  quality  partition  as  input.  With  more 
memory  available,  a high-quality  partition  was  then  calculated  using  ParMETIS.  An  al- 
ternative and  presumably  faster  approach  would  be  to  increase  the  memory  per  processor 
by  using  more  nodes.  This  was  tried  unsuccessfully  up  to  6 nodes  for  the  100  M case  (i.e. 
3 processors  per  node),  at  which  point  the  “successive  partitioning”  strategy  described 
previously  was  adopted. 

Tables  1 and  2 summarize  some  key  indices  measured  during  preprocessing.  The  most 
expensive  and  apparently  least  scalable  component  of  preprocessing  is  the  cooperative 
reordering  and  writing  of  the  partition  files  using  MPI-I/O.  In  the  worst  case,  however, 
preprocessing  of  the  100  M cv  grid  with  1024  partitions  required  4.5  hours. 

2.2.  Large-Scale  Simulation 

In  the  process  of  verification  and  validation  of  CDP-a,  a variety  of  incompressible  cold 
flow  simulations  have  been  performed.  The  largest  of  these  have  been  35  M cv  and  100  M 
cv  simulations  of  the  full  P&W  combustor  (all  18  injectors).  Figure  2 presents  a snapshot 
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Table  2.  Parallel  preprocessing  of  100  M cv  grid  - 2 cases 


exact  cv  count 
number  of  msh  files,  Nm 
poor  quality  partitions,  Np> 
final  partitions,  Np 

I/O: 

input  msh  files 

poor  quality  partition  files 

final  partition  files 

Wall  clock  time: 

1.  PPAREM  on  Nm  processors: 

read/parse  msh  files 
octree  reconnect  geometry 
poor  quality  partition 
write  poor  partition  files 

2.  PPAREM  on  Npi  processors: 

read  poor  quality  partition 
partition  (ParMETIS) 
write  Np  partition  files 
Total  wall  clock  time 


100,241,208 

100,241,208 

18 

18 

96 

96 

512 

1024 

14.28  GB 

14.28  GB 

9.65  GB 

9.65  GB 

8.86  GB 

8.78  GB 

1550  s 

1550  s 

75  s 

75  s 

2 s 

2 s 

1538  s 

1538  s 

122  s 

183  s 

387  s 

501  s 

6445  s 

12385  s 

10121  s (2.8  h) 

16236  s (4.5  h) 

Figure  2.  Geometry  and  contours  of  the  instantaneous  velocity  magnitude  on  a plane  cut 
through  the  full  P&W  6000  combustor  simulation,  35  M cv’s. 


of  this  simulation  illustrating  its  substantial  geometric  complexity  and  fidelity.  The  max- 
imum ratio  of  length  scales  (ratio  of  overall  domain  size  to  the  smallest  control  volume 
sizes)  is  about  104.  When  the  new  code  is  completed  and  contains  all  the  combustion  and 
spray  capability  of  the  existing  CDP,  the  ability  to  run  such  large  simulations  will  al- 
low the  investigation  of  important  physical  phenomena,  including  azimuthal  combustion 
instabilities.  Their  purpose  at  this  point,  however,  is  mainly  for  computer  science. 
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2.2.1.  Multigrid 

Multigrid  solvers  will  play  an  important  role  in  achieving  speed  and  scalability  for  the 
very  large  scale  simulations  made  possible  with  the  new  code.  As  part  of  the  CDP-a 
development,  an  interface  was  written  to  the  hypre  solvers  of  LLNL’s  CASC  group  (Fal- 
gout  et  al.  2002),  and  specifically  the  algebraic  multigrid  solver  (AMG).  Tests  on  large 
scale  cold-flow  simulations  (35  M cv’s  on  between  160  and  512  processors  of  Frost)  where 
the  existing  preconditioned  conjugate  gradient  solver  (PCG)  for  pressure  was  replaced 
by  AMG  yielded  a consistent  overall  speedup  of  over  4 times.  At  the  largest  problem 
size  investigated  of  100  M cv’s,  however,  it  was  not  possible  to  use  AMG  due  to  memory 
allocation  errors,  even  when  all  64  nodes  of  Frost  were  dedicated  to  this  problem  . Based 
on  these  results,  the  following  modifications  were  made  to  the  solver-related  research  in 
the  combustor  group: 

• Although  the  existing  CDP  code  will  be  phased  out  in  the  next  few  months  as  the 
combustion  and  spray  capabilities  are  implemented  and  verified  in  CDP-a,  the  CDP-to- 
hypre  interface  was  integrated  into  the  old  CDP  code,  where  it  has  resulted  in  an  overall 
speedup  of  over  2 x for  the  multiphysics  reacting  flow  simulations.  The  effect  of  AMG 
on  both  the  old  and  new  codes  is  shown  in  figure  3. 

• The  combustor  group’s  research  into  geometric  multigrid  has  been  given  lower  prior- 
ity. While  the  performance  of  geometric  multigrid  should  theoretically  be  superior  to  the 
algebraic  method,  both  methods  have  the  property  of  linear  scalability,  and  the  promising 
results  from  hypre’s  AMG  have  meant  that  combustor  group  computer  science  resources 
can  be  concentrated  on  other  non-scalable  bottlenecks  associated  with  CDP,  including 
parallel  preprocessing  and  particle/mesh  load  balancing. 

• LLNL’s  CASC  group  has  been  involved  to  help  solve  the  memory  problems  associ- 
ated with  the  very  large  scale  simulations,  and  also  to  provide  expertise  in  optimizing 
the  variety  of  AMG  settings  for  our  simulations. 

2.2.2.  Scalability 

A scalability  study  of  CDP-a  has  been  performed  on  Frost.  Figure  4 summarizes  these 
results  for  the  35  M cv  full  combustor  cold  flow  LES,  showing  acceptable  scalability  on 
up  to  480  processors. 

Although  the  full  combustor  simulation  was  also  run  with  100  M cv’s  on  up  to  64  nodes 
of  Frost,  this  was  only  possible  with  the  significantly  slower  PCG  solver  for  pressure. 
Investigations  of  the  code  scalability  at  these  larger  problem  sizes  and  over  a greater 
range  of  processors  will  be  made  once  the  memory  problems  associated  with  the  AMG 
solver  have  been  solved. 

2.2.3.  Code  Performance 

Table  3 gives  some  of  the  key  computer  science  indices  for  the  new  CDP-a  code.  The 
total  memory  requirement  and  total  file  size  scale  linearly  with  problem  size,  so  these 
indices  are  reported  per  M cv’s.  We  note  that  the  AMG  memory  estimate  is  based  on 
a number  of  smaller  simulations  made  on  32  processors  or  fewer.  Clearly  such  an  aggre- 
gate memory  requirement  does  not  give  the  complete  picture  for  large-scale  simulations. 
Assuming  linear  scaling,  the  100  M cv  simulation  should  require  100  M x 3.5  GB  = 350 
GB,  less  than  the  physical  memory  available  on  32  and  certainly  64  nodes  of  Frostf.  As 
discussed  earlier,  however,  attempts  to  run  the  100  M cv  simulation  with  AMG  on  both 
32  and  64  nodes  of  Frost  resulted  in  memory  allocation  errors. 

f Frost  has  64  nodes  with  16  GB/node  and  16  processors/node 
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H scalars 

■ momentum  + pressure 
□ sprays 


Q momentum 
■ pressure 


wall  clock  time  per  time  step  [seconds] 


wall  clock  time  per  time  step  [seconds] 


(a)  old  CDP,  multiphysics  LES  (b)  new  CDP-a,  cold  flow  LES 

FIGURE  3.  Speedup  of  the  old  and  new  CDP  codes  by  changing  the  pressure  solver  from  pre- 
conditioned conjugate  gradient  (PCG)  to  hypre’s  algebraic  multigrid  (AMG):  a)  Old  CDP  code, 
P&W  reacting  flow  single  sector  simulation  with  2 M cv’s  on  80  processors  of  Frost  («  25,000 
cv’s/processor).  b)  New  CDP-a,  cold  flow  LES  of  full  combustor,  35  M cv’s  on  480  processors 
of  Frost  (ss  73, 000  cv’s/processor). 


Processors,  Np 

FIGURE  4.  Scalability  of  the  new  CDP-a  code  on  Frost  (with  hypre  AMG  solver  for  pressure). 
Problem  is  the  cold  flow  LES  of  the  full  combustor  geometry  with  35  M cv’s.  Scaling  is  reported 
relative  to  the  160  processor  case,  which  was  the  smallest  number  of  Frost  processors  that  could 
run  this  problem. 


2.3.  Particle/Mesh  Load  Balancing 

The  geometric  locality  of  particles  combined  with  the  sequential  nature  of  particle/mesh 
time  advancement  can  create  load  imbalance.  By  integrating  components  of  the  parallel 
preprocessor  with  CDP-a,  it  was  straight-forward  to  add  particle/mesh  load  balancing 
capabilities  to  the  new  code.  Because  the  code  has  very  long  run  times  relative  to  the 
startup  time  (i.e.  reading  the  partition  and  restart  files),  we  developed  the  load  balancing 
strategy  to  work  like  a restart.  When  a threshold  value  of  load  imbalance  is  reached,  a 
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Table  3.  Some  computer  science  indices  for  CDP-q 


Total  memory  requirement 

- with  PCG  solver  for  pressure 

1.5  GB  per  M cv’s 

- with  AMG  solver  for  pressure 

3.5  GB  per  M cv’s 

I/O  (MPI-2): 

- input  partition  files 

0.10  GB  per  M cv’s 

- input  restart  or  output  result  files 

0.20  GB  per  M cv’s 

Floating  point  performance  (Frost): 

- sustained  Megaflops  per  processor 

78.9  Mflops 

- percentage  of  peak 

5.3  % 

- flops  to  memory  references 

0.294 

Sample  wall  clock  time  (35  M cv  on 

480  processors  of  Frost,  cold  flow): 

- startup  (read  partition  & restart  files) 

13  min 

- 20,000  time  steps  at  11  s/step 

61  h 

- writing  result/restart  files 

5 min 

- cooperative  writing  2D  plane  data  file 

1 min 

new  partition  is  calculated  using  the  multi-constraint  option  available  in  ParMETIS. 
This  new  partition  will  have  (approximately)  the  same  number  of  cv’s  and  particles 
on  each  processor.  This  new  partition  is  then  cooperatively  written  to  files  using  the 
routines  developed  for  the  parallel  preprocessor.  The  simulation  is  then  stopped,  and  can 
be  restarted  using  the  new  partition. 

At  the  time  of  this  writing,  these  capabilities  have  been  tested  on  small  problems  using 
model  particle  cost  distributions  only.  Figure  5 illustrates  a model  particle  cost  function 
and  an  initial  32-partitioning  produced  by  the  parallel  preprocessor  using  ParMETIS  in 
single-constraint  mode.  ParMETIS  was  then  used  to  repartition  the  problem  with  the 
dual  constraints  of  equal  numbers  of  cv’s  and  particles  on  each  partition.  The  theoretical 
simulation  timings  for  the  initial  partition  and  multi-constraint  partition  are  shown  in  fig- 
ure 6 a)  and  b)  respectively.  Times  are  normalized  relative  to  the  perfectly  balanced  case. 
For  this  particular  particle  distribution,  the  initial  partition  has  a normalized  computa- 
tion time  of  over  4.  Using  a multi-constraint  repartitioning,  the  normalized  computation 
time  is  reduced  to  1.04.  This  implies  a speedup  of  about  4 x,  assuming  no  change  in 
parallel  efficiency.  The  multi-constrained  partitioning  does,  however,  increase  the  edge 
cut  - in  this  case  by  about  40%.  The  extent  to  which  this  reduces  the  theoretical  speedup 
will  be  determined  when  this  load  balancing  capability  is  tested  on  real  problems. 


3.  Spray  Models  in  CDP 

The  objective  of  the  spray  simulation  effort  is  to  develop  a numerical  framework  based 
on  Lagrangian  particle  models  to  perform  high-fidelity  LES  of  reacting  multiphase  flows 
encountered  in  realistic  gas-turbine  combustion  chambers.  Three  regimes  of  flow  devel- 
opment are  commonly  observed  inside  these  combustors  (as  shown  in  Fig.  7):  1)  ’primary 
breakup  regime’  and  ‘dense  regime’,  where  the  liquid  film  exhibits  large  scale  coherent 
structures  that  interact  with  the  gas-phase  and  disintegrate  into  filaments  by  forming 
Kelvin-Helmholtz  type  instability,  2)  ‘intermediate  regime’  where  the  liquid  blobs  formed 
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(a)  model  particle  cost  (b)  initial  32-partitioning 


Figure  5.  Sample  problem  used  to  demonstrate  particle/mesh  load  balancing  capabilities  in 

CDP-a. 


■gas  didle  Opartlde  Oldie  Bgas  Oldie  Opartlcle  Oldie 


(a)  single-constraint  (b)  multiple-constraint  32- 

32-partition,  edge  cut  = partition,  edge  cut  = 23043, 

18038,  normalized  time  = normalized  time  = 1.04 

4.05 

Figure  6.  Computation  time  on  each  processor  for  sample  particle/mesh  load  balancing 

problem  shown  in  figure  5. 

undergo  secondary  breakup,  and  3)  ‘dilute  regime’  where  the  droplets  evaporate,  the  fuel 
vapor  mixes  with  the  surrounding  hot  gas  giving  turbulent  spray  flames.  Droplet  defor- 
mation and  collision  are  also  important  features  in  the  intermediate  regime.  Figure  7 
summarizes  the  characterization  of  the  spray  formation  into  the  three  regimes,  based  on 
the  ratio  of  the  characteristic  length  scale  of  the  liquid  to  the  available  grid  resolution 
if  Ax  and  the  liquid  phase  volume  fraction  ©p. 

The  main  task  was  to  integrate  models  capturing  these  complex  phenomena  into  the 
unstructured  LES  code  (CDP),  perform  comprehensive  validation  simulations  of  each 
model,  and  initiate  large-scale  turbulent,  multiphase,  reacting  flow  simulation  in  realistic 
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dilute  regime 


intermediate  regime 


primary  breakup 


dense  regime 
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FIGURE  7.  Regimes  of  spray  formation  defined  by  the  ratio  of  the  characteristic  length  scale  of 
the  liquid,  i,  to  the  local  grid  spacing,  Ax  and  the  liquid  phase  volume  fraction,  ©p. 

combustor  geometry.  A summary  of  Eulerian/Lagrangian  equations  with  two-way  cou- 
pling for  interphase  mass,  momentum,  and  energy  transport  is  given  below.  The  subgrid 
models  for  sprays  including  secondary  breakup,  droplet  deformation  and  drag,  droplet 
evaporation,  a hybrid  particle-parcel  approach  for  sprays,  and  accounting  for  finite-size 
effect  of  droplets  in  the  dense  spray  regime  are  described  in  brief.  Also,  two  novel  ap- 
proaches to  more  accurately  simulate  the  initial  ’primary  breakup  regime’  and  the  ’dense 
regime’  are  briefly  outlined. 

3.1.  Gas-Phase  Equations 

We  solve  the  variable  density,  low-Mach  number  equations  with  two-way  coupling  be- 
tween the  gas-phase  and  liquid  particles.  The  formulation  is  based  on  the  flamelet /progress 
variable  approach  developed  by  Pierce  & Moin  (2001)  for  LES  of  non-premixed,  turbulent 
combustion.  The  gas-phase  continuity,  scalar,  and  momentum  equations  are, 

^ = <3-’> 

<-> 


&p-gC  dJTgCui  d dC  dqCj  -r- 

~dT + -&r  - 6^ac~}  - — + 


d/V^»  Qp^UjUj  __dp_  d{2nSjj)  _ dq,3  y 

On  Ok  Ok  <"  Ok  r\  ‘-'l 


pg  = p-g(Z,C,Z”)  (3.5) 

- _ 1 dui  duj  1 duk 

Sij~2  {duj  + dui)  Z5ijd7k  (3'6) 

where  the  unclosed  transport  terms  in  the  momentum  and  scalar  equations  are  grouped 
into  the  residual  stress  qig,  and  residual  scalar  flux  qz3andqc3.  The  dynamic  Smagorinsky 
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model  by  Moin  et  al.  (1991)  is  used  to  close  these  subgrid  terms  as  demonstrated  by 
Pierce  &i  Moin  (2001).  For  a two-fluid  (air  + fuel)  mixture,  one  conserved  scalar  (the 
mixture  fraction,  Z)  and  a non-conserved  scalar  (the  progress  variable,  C)  are  solved.  The 
gas-phase  properties  such  specific  heat,  molecular  weight,  density,  viscosity,  heat-release, 
and  source  terms  in  the  progress- variable  equation,  ujc,  are  obtained  from  lookup  tables 
generated  using  flamelet  theory  for  non-premixed  combustion  and  are  dependent  on  the 
local  values  of  Z,  C,  and  the  subgrid  mixture  fraction  fluctuations,  Z" . 

With  the  presence  of  the  liquid  phase  and  inter-phase  mass,  momentum,  and  energy 
transport,  additional  source  terms  are  added  into  the  continuity,  scalar  transport  and 
momentum  equations.  As  the  droplets  evaporate  the  heat  of  vaporization  is  taken  from 
the  gas-phase  and  there  is  evaporative  cooling  of  the  surrounding  gas.  This  gives  rise  to  a 
sink  term  in  the  energy  equation.  By  assuming  adiabatic  walls  and  unity  Lewis  number, 
the  energy  and  scalar  equations  have  the  same  boundary  conditions  and  are  linearly 
dependent.  Only  one  scalar  equation  (for  mixture  fraction,  Z)  is  solved  and  other  scalars 
including  temperature  are  deduced  using  flamelet  tables.  The  evaporative  cooling  effect 
(heat  of  vaporization)  is  accounted  for  in  the  equation  of  state  during  the  generation  of 
the  flamelet  tables.  The  heat  content  of  the  liquid  fuel  is  taken  into  account  by  computing 
an  effective  gaseous  fuel  enthalpy  and  is  used  in  solving  the  flamelet  equations. 

3.2.  Modeling  Primary  Breakup/Dense  Regime 

The  initial  breakup  of  the  turbulent  liquid  film  into  large  coherent  structures  is  called 
primary  breakup.  It  is  dominated  by  the  interaction  of  the  turbulent  liquid  film  with 
the  surrounding  turbulent  gas-phase  giving  rise  to  liquid  surface  instability  waves.  These 
interfacial  instabilities  are  important  in  the  overall  spray  evolution  and  droplet  formation 
process.  However,  the  dynamics  of  the  phase  interface  is  highly  complex  and  as  of  today 
not  well  understood.  It  in  fact  remains  one  of  the  outstanding  unresolved  problems  in  the 
area  of  spray  simulation  how  to  model  this  initial  breakup  in  turbulent  environments  cor- 
rectly. The  majority  of  the  atomization  modeling  effort  is  based  on  Lagrangian  particle 
tracking  and  secondary  breakup  mechanisms  as  described  later.  Implicit  in  these  models 
are  the  assumptions  that  a)  the  characteristic  length  scale  of  the  liquid,  £,  is  small  com- 
pared to  the  grid  size,  Ax,  and  b)  the  liquid  volume  fraction  ©p  is  small.  Although  these 
models  are  successful  in  predicting  secondary  breakup  in  the  ‘intermediate’  and  ‘dilute’ 
regimes,  they  are  not  applicable  in  the  ‘primary  breakup’  and  ‘dense’  regime,  because 
their  basic  assumptions  are  invalid. 

Two  different  but  in  essence  complementary  approaches  are  being  pursued  within  the 
group  to  correctly  predict  the  ‘primary  breakup’  and  ‘dense’  regime.  On  the  one  hand, 
a novel  Eulerian  based  Large  Surface  Structure  model  is  being  developed  that  explicitly 
tracks  the  larger  scale  dynamics  of  the  phase  interface,  If  Ax  > 1.  On  the  other  hand, 
Lagrangian  spray  models  are  being  extended  to  the  ‘dense  regime’  by  accounting  for 
drop/particle  sizes  that  are  comparable  to  but  still  smaller  than  the  grid  size,  If  Ax  < 1. 

3.2.1.  Large  Surface  Structure  Model 

The  objective  of  the  Large  Surface  Structure  (LSS)  Model  is  to  correctly  capture  the 
dynamics  of  the  phase  interface  in  the  primary  breakup  regime.  To  this  end,  we  propose  to 
follow  in  principle  the  LES  approach  for  turbulent  flows.  All  interface  dynamics  occurring 
on  scales  larger  than  the  grid  size  are  explicitly  resolved  and  captured  by  a Eulerian  level 
set  approach,  whereas  interface  dynamics  occurring  on  the  subgrid  scales  are  described 
by  an  appropriate  subgrid  model.  Details  of  this  model  and  some  preliminary  test  cases 
and  results  are  given  in  a separate  paper  in  this  annual  research  brief  (Herrmann  2003). 
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3.2.2.  Modeling  Finite-Size  Droplets 

In  a parallel  effort,  the  Lagrangian  framework  is  extended  to  account  for  drop/particle 
sizes  that  are  comparable  to  but  still  smaller  than  the  grid  size,  ij  Ax  < 1.  The  equations 
are  based  on  the  original  spray  formulation  developed  by  Dukowicz  (1980)  which  consists 
of  Eulerian  fluid  and  Lagrangian  particle  calculations,  and  accounts  for  the  displacement 
of  the  fluid  by  the  particles  as  well  as  the  momentum  interchange  between  them.  The 
gas-phase  equations  given  above  are  modified  to  take  into  account  the  volume  fraction 
(0  = 1 — ©p),  where  ©p  represents  the  volume  occupied  by  the  liquid-phase  in  a given 
control  volume.  The  gas-phase  density,  pg  and  pressure  p are  then  replaced  by  (ps0)  and 
(pO)  in  the  gas-phase  governing  equations. 

In  brief,  the  Lagrangian  particles  are  advanced  using  Lagrangian  particle  tracking. 
The  particle  volume  fraction  0p  is  then  computed  on  each  grid  point  and  then  the 
gas-phase  equations  are  advanced.  Use  of  volume  fraction,  however,  adds  to  numerical 
complexity  and  accurate  numerical  schemes  conserving  liquid  mass  are  being  developed. 
Apte  et  al.  (2003c)  summarize  the  theoretical  formulation  and  present  some  preliminary 
test  cases  in  particulate  flows. 

3.3.  Modeling  Intermediate  Regime 

The  particle  motion  is  simulated  using  the  Basset-Boussinesq-Oseen  (BBO)  equations 
(Crowe  et  al.  1998).  It  is  assumed  that  the  density  of  the  particle  is  much  larger  than 
that  of  the  fluid  (pp/pg  ~ 103) , particle-size  is  small  compared  to  the  turbulence  integral 
length  scale,  and  that  the  effect  of  shear  on  particle  motion  is  negligible.  The  high  value 
of  density  ratio  implies  that  the  Basset  force  and  the  added  mass  term  are  small  and  are 
therefore  neglected.  Under  these  assumptions,  the  Lagrangian  equations  governing  the 
droplet  motions  become 


^=Drl  »-u,)  + (l-g)g  (3.8) 

where  mp  is  the  mass  of  the  droplet,  up  the  particle  velocity  components,  u gas-phase 
velocities  interpolated  to  the  particle  location,  pv  and  pg  the  particle  and  gas-phase  den- 
sities, and  g the  gravitational  acceleration.  The  drag  force  on  a solid  particle  is  modeled 
using  a drag-coefficient,  Cd, 


D. 


P solid  qpd 


Pg  I tig  up 
Pp  dp 


where  Cd  is  obtained  from  the  nonlinear  correlation  (Crowe  et  al.  1998) 


(3.9) 


cd  = Jfe(l  + aR4)-  (3-10) 

Here  Rep  = dp\ug  - up\/pg  is  the  particle  Reynolds  number.  The  above  correlation 
is  applicable  for  Rep  < 800.  The  constants  a = 0.15,6  = 0.687  yield  the  drag  within 
5%  from  the  standard  drag  curve.  Modifications  to  the  solid  particle  drag  are  applied  to 
compute  drag  on  a liquid  drop  and  are  given  below. 
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3.4.  Deformation  and  Drag  Models 

In  order  to  quantify  the  effect  of  droplet  deformation  on  drag,  Helenbrook  Sz  Edwards 
(2002)  performed  detailed  resolved  simulations  of  axisymmetric  liquid  drops  in  uniform 
gaseous  stream.  The  simulations  were  performed  using  an  hp- finite  element  method  (He- 
lenbrook 2002).  An  unstructured  mesh  of  triangles  which  deforms  with  the  interface  along 
with  a dynamic  mesh  adaptation  algorithm  was  used  allowing  higher-order  accuracy  to 
be  obtained  even  though  there  is  a discontinuity  at  the  interface.  This  gives  results  for 
the  drag  which  are  accurate  to  within  1%.  Based  on  their  computations  for  a range  of 
density  and  viscosity  ratios,  range  of  Weber  (We),  Ohnesorge  (Oh),  and  Reynolds  num- 
bers (Re),  a correlation  was  developed  that  provides  the  amount  of  deformation  in  the 
form  of  ellipticity,  E,  which  is  defined  as  the  ratio  of  the  height  to  width  of  the  drop. 
This  is  given  as, 


E = 1 - O.llWe0'82  + 0.013. — Oh~0  55 We11  (3.11) 

\ Pg  Ml 

where  pi,  pg  are  the  viscosities  of  the  liquid  and  gas  and  pi,  pg  are  the  densities,  re- 
spectively. The  non-dimensional  Weber  and  Ohnesorge  numbers  axe  defined  as,  We  = 
pgU2dp/o  and  Oh  = pi/y/piadp,  where  U is  the  relative  velocity  between  the  gas  and 
liquid,  dp  the  diameter  of  the  droplet,  and  a the  surface  tension.  Accordingly,  E < 1 
indicates  that  the  drops  have  more  width  than  height  with  deformation  in  the  direction 
perpendicular  to  the  relative  velocity.  These  shapes  are  called  oblate  shapes.  Similarly, 
E > 1 gives  elongation  in  the  direction  of  the  relative  velocity  giving  rise  to  prolate 
shapes.  E = 1 implies  spherical  shapes. 

The  effect  of  droplet  deformation  is  to  change  the  drag  force.  This  effect  is  modeled 
by  using  an  effective  equatorial  droplet  diameter,  d*  = dpE~1/3.  The  particle  Reynolds 
number  is  also  modified,  Re*  = RepE~1/>3  (Helenbrook  & Edwards  2002).  This  is  used 
in  equations  (3.9,  3.10)  to  obtain  the  modified  drag.  In  addition  the  effect  of  internal 
circulation  is  modeled  by  changing  the  drag  on  a solid  sphere  as 

3.5.  Stochastic  Model  for  Secondary  Breakup 

Performing  simulations  of  primary  atomization  where  one  tracks  the  liquid-gas  interface 
in  realistic  combustor  geometries  is  computationally  intensive.  Development  of  numeri- 
cal techniques  based  on  hybrid  Eulerian/Lagrangian  Level-set /Particle  tracking  schemes 
to  capture  the  primary  atomization  process  is  ongoing  and  will  be  implemented  into 
CDP.  However,  the  current  status  is  to  compute  the  atomization  process  using  advanced 
stochastic  secondary  breakup  models  developed  (Apte  et  al.  2003).  Emphasis  is  placed 
on  obtaining  the  correct  spray  evolution  characteristics  such  as  liquid  mass  flux,  spray 
angle,  and  droplet  size  distribution.  The  liquid  film  is  approximated  by  large  drops  with 
size  equal  to  the  nozzle  diameter  and  undergoes  deformation  and  breakup.  The  effect  of 
high  mass-loading  on  the  gas-phase  momentum  transport  is  captured  through  two-way 
coupling  between  the  two  phases. 

In  this  model,  the  characteristic  radius  of  droplets  is  assumed  to  be  a time-dependent 
stochastic  variable  with  a given  initial  size-distribution.  The  breakup  of  parent  blobs  into 
secondary  droplets  is  viewed  as  the  temporal  and  spatial  evolution  of  this  distribution 
function  around  the  parent-droplet  size  according  to  the  Fokker-Planck  (FP)  differen- 
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tial  equation.  This  distribution  function  follows  a certain  long-time  behavior,  which  is 
characterized  by  the  dominant  mechanism  of  breakup: 


VT(*,t)+mnpt)  = 1 (fJ) 


82T(x,t) 


(3.13) 


dt  ' dx  ' dx2  ' 

where  the  breakup  frequency  (v)  and  time  (t)  are  introduced.  Here,  T(x,t ) is  the  dis- 
tribution function  for  x — log(rj),  and  r}  is  the  droplet  radius.  Breakup  occurs  when 
t > tbreakup  = l/"-  The  value  of  the  breakup  frequency  and  the  critical  radius  of  breakup 
are  obtained  by  the  balance  between  the  aerodynamic  and  surface  tension  forces.  The  sec- 
ondary droplets  are  sampled  from  the  analytical  solution  of  equation  (3.13)  corresponding 
to  the  breakup  time-scale.  The  parameters  encountered  in  the  FP  equation  ((£)  and  (£2)) 
are  computed  by  relating  them  to  the  local  Weber  number  for  the  parent  blob,  thereby 
accounting  for  the  capillary  forces  and  turbulent  properties.  As  new  droplets  are  formed, 
parent  droplets  are  destroyed  and  Lagrangian  tracking  in  the  physical  space  is  continued 
till  further  breakup  events. 


3.6.  Modeling  Dilute  Regime 

Typical  spray  simulations  do  not  resolve  the  temperature  and  species  gradients  around 
each  droplet  to  compute  the  rate  of  evaporation.  Instead,  evaporation  rates  are  estimated 
based  on  quasi-steady  analysis  of  a single  isolated  drop  in  a quiescent  environment  (Faeth 
1977,  Faeth  1983).  Multiplicative  factors  are  then  applied  to  consider  the  convective  and 
internal  circulation  effects.  We  model  the  droplet  evaporation  based  on  a ‘uniform-state’ 
model.  The  Lagrangian  equations  governing  particle  mass  and  heat  transfer  processes 
are  well  summarized  by  Oefelein  (1997)  and  are  described  in  brief. 


d 

— (mp)  = -rhp  (3.14) 

mpCPi~lt  (Tp)  - hvnd2(Tg  - Tp)  - fnpAhv  (3.15) 

where  A hv  is  the  latent  heat  of  vaporization,  mp  mass  of  the  droplet,  Tp  temperature  of 
the  droplet,  and  Cpi  the  specific  heat  of  liquid.  The  diameter  of  the  droplet  is  obtained 
from  its  mass,  dp  = {&mp/,npp)1^ . hp  is  the  effective  heat-transfer  coefficient  and  is 
defined  as 


K = K (f?)  /(Tg-Ts)  (3.16) 

where  ks  is  the  effective  conductivity  of  the  surrounding  gas  at  the  droplet  surface. 
The  subscript  ‘s’  stands  for  the  surface  of  the  droplet.  The  solutions  to  the  mass  and 
temperature  equations  for  quiescent  medium  are  obtained  by  defining  Spalding  mass  and 
heat  transfer  numbers  and  making  use  of  the  Clausius-Clapeyron’s  equilibrium  vapor- 
pressure  relationship  (Faeth  1977).  In  addition,  several  convective  correction  factors  are 
applied  to  obtain  spray  evaporation  rates  at  high  Reynolds  numbers  (Faeth  1977,  Faeth 
1983). 


3.7.  Hybrid  Particle- Parcel  Technique  for  Spray  Simulations 

Performing  spray  breakup  computations  using  Lagrangian  tracking  of  each  individual 
droplet  gives  rise  to  a large  number  of  droplets  (w  10-30  million)  very  close  to  the  injec- 
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tor.  For  LES,  simulating  all  droplet  trajectories  would  be  ideal,  however,  this  gives  rise 
to  severe  load-balancing  problems  on  parallel  processors.  We  have  developed  a hybrid 
particle-parcel  scheme  to  effectively  reduce  the  number  of  particles  tracked.  A parcel  or 
computational  particle  represents  a group  of  droplets,  Npar,  with  similar  characteristics 
(diameter,  velocity,  temperature).  During  injection,  new  particles  added  to  the  compu- 
tational domain  are  pure  drops  ( Npar  = 1).  These  drops  move  downstream  and  undergo 
breakup  according  to  the  secondary  breakup  model  and  produce  new  droplets. 

The  basic  idea  behind  the  hybrid-approach,  is  to  collect  all  droplets  in  a particular 
control  volume  and  group  them  into  bins  corresponding  to  their  size  and  other  properties 
such  as  velocity,  temperature  etc.  The  droplets  in  bins  are  then  used  to  form  a parcel 
by  conserving  mass,  momentum  and  energy.  The  properties  of  the  parcel  are  obtained 
by  mass-weighted  averaging  from  individual  droplets  in  the  bin.  For  this  procedure,  only 
those  control  volumes  are  considered  for  which  the  number  of  droplets  increases  above 
a certain  threshold  value.  The  number  of  parcels  created  would  depend  on  the  number 
of  bins  and  the  threshold  value  used  to  sample  them.  The  parcel  thus  created  then 
undergoes  breakup  according  to  the  above  stochastic  sub-grid  model,  however,  does  not 
create  new  parcels.  On  the  other  hand,  Npar  is  increased  and  the  diameter  is  decreased  by 
mass-conservation.  This  strategy  reduces  the  total  number  of  computational  particles  in 
the  domain.  In  a real  simulation  with  breakup  and  evaporation  models,  all  the  particles 
are  clustered  in  a very  small  region  close  to  the  injector.  This  may  still  give  rise  to  load- 
imbalance  as  only  a few  processors  would  solve  the  Lagrangian  particle  equations.  Domain 
decomposition  methods  taking  into  account  this  load  imbalance  are  being  developed  and 
will  be  tested  in  the  present  large-scale  multi-physics  simulation. 

4.  Validation  Studies  using  CDP 

A systematic  and  extensive  validation  effort  of  CDP  code  and  the  Lagrangian  spray 
modules  was  initiated  earlier  and  several  simulations  were  performed  to  compare  the 
LES  predictions  of  gas  and  liquid/solid  phases  with  the  available  experimental  data. 
Apte  et  al.  (2003a)  first  performed  an  LES  of  particle-laden  swirling,  co-annular  jet  and 
compared  the  results  with  the  experiments  of  Sommerfeld  k.  Qiu  (1991).  This  validated 
the  Lagrangian  particle  tracking  scheme  as  well  as  the  accuracy  of  the  numerical  method 
in  swirling  flows.  In  addition,  development  & validation  of  the  secondary  breakup  model 
in  simplified  combustor  geometries  was  also  completed  (Apte  et  al.  2003).  This  spray 
model  validation  effort  was  continued  and  applied  to  study  model  predictions  for  spray 
evaporation  and  breakup  in  complex  geometries. 

4.1.  Validation  of  Spray  Breakup  Model  in  PW  Frontend  Geometry 

The  stochastic  model  along  with  the  hybrid  particle-parcel  approach  were  used  to  sim- 
ulate spray  evolution  from  the  Pratt  and  Whitney  injector  nozzle.  Figure  8 shows  the 
instantaneous  snapshot  of  the  spray  field  along  with  the  injector  geometry.  The  experi- 
mental data  set  was  obtained  by  mounting  the  injector  in  a cylindrical  plenum  through 
which  gas  with  prescribed  mass-flow  rate  was  injected.  The  gas  goes  through  the  main 
and  guide  swirler  to  create  a swirling  jet  into  the  atmosphere.  Liquid  film  is  injected 
through  the  filmer  surface  which  forms  an  annular  ring.  The  liquid  mass-flow  rate  corre- 
sponds to  certain  operating  conditions  of  the  gas-turbine  engine.  For  this  case,  3.2M  grid 
points  were  used  with  high  resolution  near  the  injector.  The  domain  decomposition  is 
based  on  the  optimal  performance  of  the  Eulerian  gas-phase  solver  on  96  processors.  Due 
to  breakup,  a large  number  of  droplets  are  created  in  the  vicinity  of  the  injector.  With 
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Figure  8.  Spray  evolution  from  a realistic  gas-turbine  injector:  scatter  plot  shows  the  instan- 
taneous droplet  locations  in  z = 0 plane.  Large  size  droplets  are  injected  from  the  nozzle  wall 
in  an  annular  ring  to  form  a hollow  cone  spray. 


Figure  9.  Comparison  of  radial  variation  of  liquid  axial  mass  flux  at  various  axial  locations, 
P&W  RANS  with  TAB  model  , present  LES  with  stochastic  model,  o — o exper- 
imental error  bar. 

the  hybrid  approach,  the  total  number  of  computational  particles  tracked  at  station- 
ary state  is  around  3.5M,  which  represents  approximately  13M  droplets.  This  includes 
around  150,000  parcels.  The  load-imbalance  due  to  atomization  was  found  to  be  sig- 
nificant as  only  1/Zrd  of  the  processors  had  more  than  10,000  computational  particles. 
We  are  looking  into  advanced  load-balancing  techniques  to  reduce  this  computational 
overhead  due  to  sprays.  CDP-(o:)  has  the  capability  of  dynamic  load-balancing  based  on 
particle  imbalance  as  shown  earlier.  This  will  be  used  to  perform  these  simulations  to 
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Figure  10.  Instantaneous  contours  of  ISO-propyl  mass-fraction  superimposed  by  droplet 

locations  in  z = 0 plane. 

significantly  improve  the  efficiency  and  speedup.  It  will  also  facilitate  us  to  investigate 
advanced  collision/coalescence  models. 

Figure  (9)  shows  comparisons  with  the  experimental  data  of  radial  variation  of  liq- 
uid mass-flowrates  using  LES  with  stochastic  model  for  secondary  breakup.  Also  shown 
are  the  predictions  made  by  Pratt  k Whitney’s  RANS-based  simulation  with  a Taylor- 
Analogy  Breakup  (TAB)  model.  The  LES  results  are  generally  in  good  agreement  with 
the  experiments.  The  liquid  mass-flowrates  basically  relate  to  the  spray  angle  and  droplet 
size  distributions  obtained  after  breakup.  The  Sauter  Mean  Diameter  (SMD)  are  well 
captured  by  the  present  LES  simulation  (within  7%  of  the  experimental  data).  The  di- 
ameter distributions  when  compared  with  the  experimental  data  indicates  a broad  size- 
distribution  as  opposed  to  the  one  predicted  by  RANS  with  the  TAB  model.  The  size 
distribution  also  indicates  presence  of  large  number  of  small  size  droplets  which  could 
be  attributed  to  the  lack  of  models  addressing  droplet  collision  and  coalescence.  In  ad- 
dition, the  initial  droplet  size  at  the  injector  nozzle  is  assumed  to  be  constant  in  these 
simulations,  whereas  it  may  vary  depending  on  the  local  conditions  governing  primary 
atomization.  A further  investigation  with  inclusion  of  collision  models  as  well  as  using  a 
size  distribution  at  the  inlet  should  be  performed  in  order  to  accurately  predict  the  spray 
characteristics. 


4.2.  Validation  of  Spray  Evaporation  Model 

In  order  to  validate  the  evaporation  model  and  the  variable  density  formulation  in  CDP, 
simulation  of  a coaxial  non-swirling  jet  was  performed  following  the  configuration  used 
in  the  experiments  by  Sommerfeld  k Qiu  (1998).  This  flow  configuration  was  chosen 
because  of  its  direct  relevance  in  gas-turbine  combustion  chambers.  In  addition,  in  these 
experiments  the  boundary  conditions  for  the  liquid  phase  specifying  the  inlet  droplet  size 
distribution  and  their  correlation  with  droplet  velocity  is  well-defined.  The  gas-phase 
temperatures  are  not  high  enough  to  produce  spray  flames.  This  isolates  the  droplet 
evaporation  problem  from  spray  breakup  and  combustion  and  is  very  useful  in  validating 
the  evaporation  models  used  in  CDP.  Accordingly,  only  one  scalar  equation  (mixture 
fraction)  is  solved  and  the  gas-phase  temperature  is  obtained  from  the  ideal  gas  law. 

The  selection  of  droplet  size  and  velocity  at  the  injection  surface  is  based  on  the  given 
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Figure  11.  Comparison  of  radial  variation  of  evaporation  statistics  at  various  axial  locations 

with  the  experimental  data  of  Sommerfeld  & Qiu  (1998), LES,  o experimental  data:  a) 

mean  and  rms  of  axial  particle  velocity,  b)  mean  axial  liquid  mass  flux,  c)  mean  and  rms  of 
particle  droplet  diameter. 


experimental  PDF  of  size- velocity  correlations.  The  grid  used  consists  of  1.5M  hexahedral 
cells  and  around  0.5M  particles  were  present  in  the  computational  domain  at  stationary 
state.  Figure  (10)  shows  a snapshot  of  mixture  fraction  contours  superimposed  by  scatter 
plot  of  ISO-propyl  alcohol  particles  in  a coaxial  combustor.  The  droplets  are  injected 
near  the  inlet  circular  wall  of  cross-sectional  diameter  40  mm.  The  droplet  velocity- 
size  correlation  depicts  a conical  spray  with  a spray  angle  of  around  60°.  Figure  (11) 
shows  the  droplet  statistics.  The  radial  variations  of  mean  and  rms  droplet  velocity, 
mean  axial  mass  flux,  and  mean  and  rms  droplet  diameter  compared  with  experimental 
data  of  Sommerfeld  & Qiu  (1998)  are  shown.  The  profiles  are  in  good  agreement  with 
the  experimental  data.  The  droplet  mean  diameter  profile  is  typical  of  the  hollow  cone 
atomizer,  where  smaller  droplets  are  found  in  the  core  region  and  larger  droplets  near  the 
edge  of  the  spray.  Away  from  the  inlet  section,  the  droplet  mean  diameter  distribution 
becomes  more  uniform  over  the  cross-section  and  slowly  decreases  in  the  downstream 
direction  because  of  evaporation.  The  axial  mass-flux  also  decreases  towards  the  exit  due 
to  evaporation  and  is  well  captured  by  the  present  simulation.  At  x = 0 and  x = 0.786 
the  profiles  of  droplet  mass  flux  show  two-peaks  associated  to  hollow-cone  spray.  Due  to 
the  recirculation  region  downstream  of  the  center-body,  the  droplet  mass-flux  becomes 
negative.  Further  downstream  spreading  of  the  spray  is  hindered  by  the  annular  air-jets 
and  maximum  of  the  mass  flux  moves  towards  the  centerline.  These  features  are  well 
captured  by  the  present  LES  computation. 
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5.  Large  Scale  Multi-physics  Simulation  using  CDP 

A multi-scale,  multi-physics  simulation  of  turbulent  reacting  flow  in  a realistic  Pratt 
& Whitney  combustor  was  initiated.  This  simulation  includes  all  the  complex  models 
for  spray  breakup,  evaporation,  and  turbulent  combustion  described  in  section(3).  This 
simulation  will  serve  as  the  first  validation  study  of  the  reacting  multiphase  flow  solver 
(CDP)  in  complex  combustor  geometry.  Figure  (12)  shows  the  combustor  geometry  in  the 
symmetry  plane  z — 0.  The  simulation  is  performed  for  a single  injector  which  represents 
one  sector  of  the  full  combustor  containing  18  injectors.  Liquid  fuel  (Jet-A)  enters  the 
combustion  chamber  through  an  annular  ring  at  the  injector  exit.  This  liquid  film  is 
approximated  by  large  drops  of  the  size  of  the  annulus  radius.  These  drops  axe  convected 
by  the  surrounding  hot  air,  they  break,  evaporate,  and  the  fuel  vapor  thus  formed  mixes 
with  the  surrounding  air  giving  a non-premixed  spray  flame.  The  computational  grid 
consists  of  1.9M  hybrid  elements  (hexes,  pyramids,  and  tets)  with  fine  resolution  close 
to  the  injector. 

The  flamelet  library  for  Jet-A  fuel  at  gas-turbine  engine  operating  conditions,  was 
generated  by  using  a surrogate  fuel  of  80%  n-Decane  and  20%  1-2-4  tri-methyl-benzene 
which  closely  follows  the  chemical  kinetics  and  reaction  rates  of  the  Jet-A  fuel.  Around 
1000  elementary  reactions  among  100  chemical  species  were  used  to  generate  these  tables. 
The  chemical  kinetics  of  surrogate  fuel  was  compared  with  original  fuel  in  terms  of 
prediction  of  pollutants  in  laminar  flames  and  showed  good  agreement.  In  this  multiphase 
simulation,  the  flamelet  tables  are  generated  by  considering  the  heat  of  vaporization  of 
the  liquid  fuel.  This  is  obtained  by  using  heat-content  for  an  equivalent  gaseous  fuel 
which  gives  large  densities  in  the  pure  fuel-vapor  region  (Z  — 1). 

Figure  (12)  shows  instantaneous  snapshots  of  progress  variable  (C),  mixture  fraction 
( Z ),  normalized  temperature  (T/To),  and  velocity  magnitude  (urnag)  in  the  z = 0 sym- 
metry plane.  Scatter  plot  of  droplets  forming  a conical  spray  is  also  shown.  A highly 
complex  unsteady  flame  is  observed  in  this  simulation.  High  temperatures  in  the  com- 
bustion chamber  are  reduced  by  large  mass-influx  through  the  dilution  holes.  This  reduces 
the  exit  temperature  considerably.  Strong  recirculation  zone  with  conical  spray  flame  is 
observed  near  the  nozzle.  The  droplets  evaporate  completely  close  to  injector  and  do  not 
appear  beyond  ( x = 2).  The  experimental  data  available  for  validation  includes  pres- 
sure drops  across  different  components,  mass-splits  through  the  inner  and  outer  dilution 
holes,  and  swirlers.  The  comparison  of  these  quantities  with  the  experimental  data  was 
shown  to  be  within  10%  for  the  cold  flow  simulation.  Similar  results  are  obtained  for 
the  reacting  flow  case.  In  addition,  the  exit  plane  temperature  was  measured  at  5 loca- 
tions. The  temperature  field  is  within  10-15%  of  the  experimental  data  near  the  center. 
The  predicted  temperature  near  the  walls  is  much  lower  and  is  still  evolving  with  time 
indicating  that  the  flow  has  not  reached  the  statistically  stationary  state. 

This  computation  was  performed  on  80  processors  on  Frost  and  involved  around  0.5M 
computational  particles.  The  algebraic  multigrid  solver  (AMG)  developed  at  Lawrence 
Livermore  is  used  to  solve  the  Poisson  equation.  The  convergence  rate  of  the  AMG  solver 
was  compared  with  the  conjugate  gradient  solver  (PCG).  For  each  time-step,  the  conju- 
gate gradient  solver  typically  requires  2000  iterations,  whereas  the  AMG  solver  requires 
10-12  cycles  corresponding  to  500-600  conjugate  gradient  iterations.  The  AMG-solver 
gave  an  overall  speed-up  of  2-3  times  over  the  PCG  solver  and  reduced  the  computa- 
tional time  substantially. 
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6.  Summary 

CDP  is  the  flagship  LES  code  developed  by  the  combustor  group  to  perform  LES 
of  reacting  multiphase  flow  in  complex  geometry.  CDP  is  named  after  the  late  Charles 
David  Pierce  (1969-2002)  who  made  several  lasting  contributions  to  the  LES  of  reacting 
flows  and  to  the  DOE’s  ASCI  program.  CDP  is  a parallel  unstructured  finite- volume  code 
developed  specifically  for  LES  of  variable  density  low  Mach-number  flows.  It  is  written  in 
Fortran  90,  uses  MPI,  and  has  integrated  combustion  and  spray  modules.  A wide  range 
of  simulations  of  multiphase,  reacting  turbulent  flows  have  been  performed  using  CDP.  It 
has  been  shown  that  the  predictive  capability  of  CDP  for  complex  flows  in  realistic  gas- 
turbine  configurations  is  good.  A major  redesign  and  rewrite  of  CDP  has  been  initiated 
to  improve  and  add  to  the  current  capabilities.  In  the  past  year,  major  accomplishments 
included  the  following  milestones: 

• Multi-physics  simulation  in  complex  geometry:  The  first  multi-physics  simulation  in- 
cluding fuel  spray  breakup,  coalescence,  evaporation,  and  combustion  is  being  performed 
in  a single  periodic  sector  (1/184^)  of  an  actual  Pratt  & Whitney  combustor  geometry. 

• Large  scale  simulation:  Performed  large  eddy  simulations  of  the  complete  combustor 
geometry  (all  18  injectors)  with  over  100  million  control  volumes  using  CDP-a. 

Further  development  of  CDP-a  will  include  the  implementation  of  advanced  com- 
bustion and  primary  atomization  models,  and  modification  of  the  numerical  algorithm 
to  capture  acoustic  waves  and  combustion  instabilities  based  on  the  compressible  for- 
mulation of  Wall  et  al.  (2002).  Several  large-scale,  multiphysics  simulations  in  the  PW 
combustor  using  one  or  more  sectors  will  be  performed  to  assess  the  accuracy  of  the 
overall  formulation. 
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FIGURE  12.  Instantaneous  snapshots  (from  top  to  bottom)  of  progress  variable,  C,  mixture 
fraction,  Z , temperature  normalized  by  inlet  air  temperature,  T/To,  and  velocity  magnitude, 
y'u2  + v2  + w2  superimposed  with  droplet  scatter  plot  in  Pratt  & Whitney  combustor,  z = 0 
plane.  The  geometry  has  been  distorted. 


